#include <omp.h>
#include <iostream>
#include <fftw3.h>

using namespace std;

void sub(double* b, double* a, int i, fftw_plan plan){
  double *c;
  fftw_complex *tmp1, *tmp2;
  tmp1 = (fftw_complex*)malloc(sizeof(fftw_complex)*8);
  tmp2 = (fftw_complex*)malloc(sizeof(fftw_complex)*8);
 
  fftw_execute(plan);
  free(tmp1);
  free(tmp2);
  b[i] = 2.0 * a[i];
}

int main(int argc, char** argv){
  int nthreads, tid;
  int N = 10;
  double a[N], b[N];
  int chunk = 1; 
  fftw_complex *tmp1, *tmp2;
  tmp1 = (fftw_complex*)malloc(sizeof(fftw_complex)*8);
  tmp2 = (fftw_complex*)malloc(sizeof(fftw_complex)*8);
  fftw_plan plan;
  plan = fftw_plan_dft_3d(8, 1, 1,\
			  tmp1, tmp2, FFTW_BACKWARD, FFTW_ESTIMATE);
  free(tmp1);
  free(tmp2);
  
  for(int i = 0; i < N; i++)
    a[i] = (double)i;

#pragma omp parallel \
  shared(N,nthreads, a, b, chunk, plan) \
  private(tid)
  {
    tid = omp_get_thread_num();
    cerr << "Thread " << tid << " starting.. " << endl;
#pragma omp for schedule(static)
    for(int i = 0; i < N; i++){
      sub(b,a,i, plan);
      cerr << "Thread " << tid << " calculating b[" << i << "]" << endl;
    }

  }
  
  fftw_destroy_plan(plan); 
  for(int i = 0; i < N; i++){
    printf("b[%d] = %f\n", i, b[i]);
  }

}
